Channeling is a distinct class of dissolution in complex porous media

The traditional model of solid dissolution in porous media consists of three dissolution regimes (uniform, compact, wormhole)—or patterns—that are established depending on the relative dominance of reaction rate, flow, and diffusion. In this work, we investigate the evolution of pore structure using numerical simulations during acid injection on two models of increasing complexity. We investigate the boundaries between dissolution regimes and characterize the existence of a fourth dissolution regime called channeling, where initially fast flow pathways are preferentially widened by dissolution. Channeling occurs in cases where the distribution in pore throat size results in orders of magnitude differences in flow rate for different flow pathways. This focusing of dissolution along only dominant flow paths induces an immediate, large change in permeability with a comparatively small change in porosity, resulting in a porosity–permeability relationship unlike any that has been previously seen. This work suggests that the traditional conceptual model of dissolution regimes must be updated to incorporate the channeling regime for reliable forecasting of dissolution in applications like geothermal energy production and geologic carbon storage.


Geometry Analysis with Image Analysis
The grains, pores, and pore throats were extracted from each time step in the simulations using a watershed segmentation algorithm and the Euclidean distance map of the grain and pore spaces was used to identify individual  grains and pores with the boundaries between pores as throats. An example of this method with each initial geometry is shown in Fig 2. The statistics of the grain, pore, and pore throat size distributions along with the characteristic length and velocities (at Pe=1) are given in Table 1.

Geometry Analysis with Autocorrelation
Here we compute the autocorrelation of the grains and velocities for both Model A and Model B (Fig 3). Both models have an autocorrelation function that steeply decreases towards zero with lag, over a length scale equal to the grain spacing. Model A is statistically anisotropic, with an autocorrelation function with square symmetry and prominent sidelobes reflecting the underlying grid. Model B is statistically isotropic, with no sidelobes.
The autocorrelation function of the along-flow component of the velocity field is statistically anisotropic, with rectangular symmetry. The scale length in the along-flow direction typically is similar to the grain spacing but is larger (by a factor of about five) in the cross-flow direction, as is expected for channels. For Model A, the autocorrelation has sidelobes reflecting the underlying periodicity of the medium, with wavelength equal to the grain spacing. The autocorrelation for Model B is similar, but without the sidelobes.

Velocity Distributions
Here we show the velocity rendering of each simulation at porosity=0.57 for Model A and Model B (Fig 5) and the corresponding PDFs of velocity (Fig.4).
In both models we see a clear difference between the channeling and wormhole regimes, with the wormhole regime showing more of a fat tail as velocities not inside the wormhole are proportionally slower, while in the channeling regime the PDF has a higher more distinct peak more similar to the uniform regime as is expected with more preferential flow pathways throughout the model. As structural complexity increases the distribution widens for all regimes, as is expected for more diverse pore throat size distributions in more complex models.  The equations are solved using finite volume discretization over an unstructured hybrid mesh. To build the mesh, the solid surface is described using an stl image. First, a Cartesian mesh of resolution h is generated. The mesh is then snapped onto the solid surface using the snappyHexMesh utility [1], i.e. cell containing solid are then removed and replaced by hexahedral or tetrahedral cells that match the solid boundaries. An additional layer of cells of the same resolution h is then added around the solid boundary to improve the representation of the solid surface. To decide the resolution used for the initial mesh, a convergence study on porosity and permeability was conducted for Model B ( Table 2). We observe that a resolution of 3 µm offers a good compromise between accuracy and size of computational mesh. Fig. 6 shows Model B with a zoom into a pore to observe the mesh at resolution 3 µm.  Arbitrary Lagrangian Eulerian method The equations are solved using the Arbitrary Lagrangian Eulerian (ALE) method [2], implemented in GeoChemFoam (www.github.com/geochemfoam) and the full solution procedure is presented in Fig. 8. For each time-step, the mesh points are moved with velocity w, which satisfies the Laplace equations with boundary condition (Equ. (??)) ∇ · D m ∇w j = 0 j=x,y,z (1) where D m is the diffusivity of the mesh motion, w j is the j-directional component and e j is the j-directional standard basis vector. With these equations, the mesh points will track the fluid-solid interface, and the mesh motion is diffused to avoid large volume ratio between neighbor cells. However, as the mesh points are displaced, the skewness of the mesh can increase and lead to failure of the transport solver. To avoid this, the mesh's skewness is checked at the end of each time-step, and the domain is fully remeshed upon failure. After remeshing, the velocity, pressure and concentration fields are mapped to the new mesh. In addition, topological errors can appear when two faces of the same mineral grain overlap, leading to failure of the flow or transport solver. To avoid this, the faces which are fully located in a topological error are eliminated before remeshing. These collapsing faces are identified by the following condition: a face defined as faceI collapsed if a ray leading from its center following its normal vector pointing toward the solid phase meets another face defined as faceJ at a distance lower than the grid size, and faceI and faceJ do not intersect. Following this remeshing algorithm, our numerical simulations are stable and topological errors are eliminated.

Time-stepping strategy
The simulations are performed using an adaptive time-stepping strategy based on the mesh Courant-Friedrich-Lewy (CFL) number defined as where ∆t is the time-step and h is the mesh resolution. The simulation are performed using a maximum mCF L number of 0.005, which offers a good compromise between accuracy, robustness and efficiency. Fig. 7 shows a comparison of permeability evolution as a function of porosity for Model B at P e = 1, Ki = 1 between mCFL=0.005 and mCFL=0.0025.
Appendix C: Robustness of ϕ, K and L for stochastically generated micromodel The study presented in the paper is limited to one instance of each of two stochastic models (Model A and B). Future work will focus on extending the findings to any generated geometry and in particular on linking the dissolution regimes to the parameters of the stochastic distribution. For this, it would be essential that the geometrical parameters that are used in the calculation of P e and Ki, i.e. the porosity ϕ, and the pore-scale length L, vary over a range much less than an order of magnitude, so that the calculation of P e and ki are robust over different instance of the same stochastic distribution. Table 3 shows the variation of porosity and pore-scale length for 12 instances of each stochastic distribution (Model A and B). For model A, ϕ varies between 0.430 and 0.445 and L varies between 1.04 and 1.15 ×10 −4 m; for model B, ϕ varies between 0.451 and 0.473 and L varies between 1.14 and 1.41 ×10 −4 m. This shows that the calculation of P e and Ki will be robust, as ϕ and L varies on a scale much smaller than an order of magnitude.

Appendix D: Time Sequence Videos of Dissolution
Movies S1-8 show the dissolution time series for select simulations A1-A4 and B1-B4. Movie S1: Visualisation of Model A P e=0.01 Ki=0.1 evolution of porosity and concentration. The grains are gray, with the concentration field in color. The pore throats are extracted by a watershed algorithm on the Euclidean distance map of the pore space and superimposed in white. This is an example of the compact dissolution regime.
Movie S2: Visualisation of Model A P e=1 Ki=1 evolution of porosity and concentration. The grains are gray, with the concentration field in color. The pore throats are extracted by a watershed algorithm on the Euclidean distance map of the pore space and superimposed in white. This is an example of the wormhole formation dissolution regime.
Movie S3: Visualisation of Model A P e=100 Ki=10 evolution of porosity and concentration. The grains are gray, with the concentration field in color. The pore throats are extracted by a watershed algorithm on the Euclidean distance map of the pore space and superimposed in white. This is an example of the channeling dissolution regime.
Movie S4: Visualisation of Model A P e=100 Ki=0.1 evolution of porosity and concentration. The grains are gray, with the concentration field in color. The pore throats are extracted by a watershed algorithm on the Euclidean distance map of the pore space and superimposed in white. This is an example of the uniform dissolution regime.
Movie S5: Visualisation of Model B P e=0.01 Ki=0.1 evolution of porosity and concentration. The grains are gray, with the concentration field in color. The pore throats are extracted by a watershed algorithm on the Euclidean distance map of the pore space and superimposed in white. This is an example of the compact dissolution regime.
Movie S6: Visualisation of Model B P e=1 Ki=1 evolution of porosity and concentration. The grains are gray, with the concentration field in color. The pore throats are extracted by a watershed algorithm on the Euclidean distance map of the pore space and superimposed in white. This is an example of the wormhole formation dissolution regime.
Movie S7: Visualisation of Model B P e=100 Ki=10 evolution of porosity and concentration. The grains are gray, with the concentration field in color. The pore throats are extracted by a watershed algorithm on the Euclidean distance map of the pore space and superimposed in white. This is an example of the channeling dissolution regime.
Movie S8: Visualisation of Model B P e=100 Ki=0.1 evolution of porosity and concentration. The grains are gray, with the concentration field in color. The pore throats are extracted by a watershed algorithm on the Euclidean distance map of the pore space and superimposed in white. This is an example of the uniform dissolution regime.